home *** CD-ROM | disk | FTP | other *** search
/ Aminet 33 / Aminet 33 - October 1999.iso / Aminet / misc / math / TCalcStats2c.lha / TCalcStats2c / AREXX / ANOVA_TwoWay.rexx < prev    next >
Encoding:
OS/2 REXX Batch file  |  1998-09-22  |  15.8 KB  |  837 lines
  1. /* ANOVA - Two Way */
  2.  
  3. options results
  4. if ~show('P','TCALC') then do
  5.     address command 'run turbocalc:turbocalc'
  6.     address command 'waitforport TCALC'
  7.     loadflag=1
  8. end
  9. address 'TCALC'
  10. 'DEFPUBSCREEN()'
  11. /* Add-in Rexx Math Library needed for some routines */
  12. signal on syntax
  13. if ~show('l','rexxmathlib.library') then
  14.    call addlib('rexxmathlib.library',0,-30)
  15. if ~show('l','rexxreqtools.library') then
  16.    call addlib('rexxreqtools.library',0,-30)
  17. if ~show('l','rexxsupport.library') then
  18.    call addlib('rexxsupport.library',0,-30)
  19.   /* add to library list */
  20. signal off syntax
  21.  
  22. /* Start Main Routine */
  23. if loadflag=1 then 'Load()'
  24. 'ActivateWindow()'
  25. range=rtgetstring(,"Enter Cell Range for Input","Input Request") /*,,'rt_pubscrname="TCALC"')*/
  26. colon=pos(":",range)
  27. if colon=0 then do
  28.     'Message "Please select a range before executing this script"'
  29.     'DEFPUBSCREEN("Workbench")'
  30.     exit
  31. end
  32.  
  33. /* Find cell references and cell, column numbers */
  34. start_cell=substr(range,1,colon-1)
  35. end_cell=substr(range,colon+1)
  36. start_row=cellrow(start_cell)
  37. end_row=cellrow(end_cell)
  38. start_col=cellcol(start_cell)
  39. end_col=cellcol(end_cell)
  40. NRows=end_row-start_row+1
  41. NCols=end_col-start_col+1
  42.  
  43. /* Get Number of Rows per Sample */
  44. rows=0
  45. rows=rtgetlong(,"Enter Number of Rows in each sample","Get Rows") /*,,'rt_pubscrname="TCALC"')*/
  46. if rows=0 then exit
  47. if datatype(left(rows,1),'n')=0 then do
  48.     'Message "Invalid row number"'
  49.     'DEFPUBSCREEN("Workbench")'
  50.     exit
  51. end
  52.  
  53. /* Get cell reference for output range */
  54. out_cell=rtgetstring(,"Enter Cell Reference for Output","Input Request") /*,,'rt_pubscrname="TCALC"')*/
  55. if out_cell="" then do
  56.     'DEFPUBSCREEN("Workbench")'
  57.     exit
  58. end
  59. if length(out_cell)<2 | datatype(left(out_cell,1),'n')=1 then do
  60.     'Message "Invalid cell reference"'
  61.     'DEFPUBSCREEN("Workbench")'
  62.     exit
  63. end
  64. /* Suppress Screen Redraw to Speed Things Up */
  65. 'Refresh 0'
  66.  
  67. /* Open a small output window on tcalc screen*/
  68. fo=0
  69. CR='0a'x
  70. DisplayMsg="Calculating...Please Wait."||CR||"User input is disabled during calculation."||CR
  71. if open(6Info, 'con:100/0/450/80/Progress/SCREEN TCALC', w) then do
  72.      call writeln(6Info, DisplayMsg)
  73.     fo=1
  74. end
  75. else do
  76.     'Message "TCALC Screen not available for Progress messages"'
  77. end
  78. CALL DELAY(150)
  79.  
  80. /* Get cell references for top cell in each column */
  81. 'SelectCell' start_cell
  82. do col=start_col to end_col
  83.     'GetCursorPos'
  84.     top_cell.col=result
  85.     'Column 1'
  86. end
  87. /* Get labels for later use on output */
  88. 'SelectCell' start_cell
  89. 'Column 1' 
  90. 'GetValue'
  91. testlabel=result
  92. testlabel=strip(testlabel)
  93. if datatype(testlabel,'n')=1 then do
  94.     labelflag=0
  95.     do x=start_col+1 to end_col
  96.     title.x="Column "||x
  97.     end
  98. end
  99. else do
  100.     labelflag=1
  101.     NRows=NRows-1
  102.     do x=start_col+1 to end_col
  103.     'GetValue'
  104.     str=result
  105.     title.x=translate(strip(str),"_"," ")
  106.     'Column 1'
  107.     end
  108. end
  109. 'SelectCell' start_cell
  110. w=0
  111. NR=NRows
  112. if labelflag=1 then NRows=Nrows+1
  113. do x=1 to NRows
  114.     'GetValue'
  115.     testlabel=result
  116.     testlabel=strip(testlabel)
  117.     if datatype(testlabel,'a')=1 then do
  118.         w=w+1
  119.         'GetValue'
  120.         sampletitle.w=result
  121.         say sampletitle.w
  122.     end
  123.     'CursorDown 1'
  124. end
  125. NRows=NR
  126. NumSamples=w /* Number of samples */
  127. if fo then call writech(6Info,"Progress...10 ")
  128.  
  129. /* Get data from cell range */
  130. col=start_col
  131. s=0
  132. z=0
  133. lav=0
  134. tot=0
  135. data.=0
  136. count.=0
  137. total.=0
  138. cell.=0
  139. do x=2 to NCols
  140.         s=1
  141.         col=col+1
  142.         'SelectCell' top_cell.col
  143.         if labelflag=1 then 'CursorDown 1'
  144.         do y=1 to NRows
  145.             z=z+1
  146.             'GetValue'
  147.             valtest=result
  148.             if datatype(valtest)='NUM' then do
  149.                 'GetValue'
  150.                 val=result
  151.                 data.s.x.y=val
  152.                 tot=tot+val
  153.                 total.x=tot
  154.                 count.x=1+count.x
  155.                 cell.s.x=(cell.s.x)+val
  156.             end
  157.             'CursorDown 1'
  158.             if z=rows then do
  159.                 z=0
  160.                 s=s+1
  161.             end
  162.         end
  163.     tot=0
  164.     lav=0
  165.     val=0
  166. end
  167. rowtot.=0
  168. do s=1 to NumSamples
  169.     do x=2 to NCols
  170.         rowtot.s=(cell.s.x)+rowtot.s
  171.     end
  172. end
  173. scount=rows*(NCols-1)
  174. if fo then call writech(6Info,"20 ")
  175.  
  176. /* Calculate Means */
  177. meancol.=0
  178. meanrow.=0
  179. meanrc.=0 /* Means within cells of sample */
  180. GM=0 /* Grand Mean of all observations */
  181. NumObs=0 /* Total number of observations */
  182. GT=0 /* Grand Total of all observations */
  183. do x=2 to NCols
  184.     meancol.x=total.x/count.x
  185.     NumObs=NumObs+count.x
  186.     GT=GT+total.x
  187. end
  188. do s=1 to NumSamples
  189.     meanrow.s=(rowtot.s)/scount
  190. end
  191. GM=GT/NumObs
  192. do s=1 to NumSamples
  193.     do x=2 to NCols
  194.         meanrc.s.x=(cell.s.x)/rows
  195.     end
  196. end
  197. if fo then call writech(6Info,"30 ")
  198.  
  199. /* Calculate Standard deviation and Variance */
  200. dat=0
  201. meenx=0
  202. sum.=0
  203. sd.=0
  204. var.=0
  205. z=0
  206. N=rows-1
  207. do x=2 to NCols
  208.         s=1
  209.         do y =1 to NRows
  210.             meen=meanrc.s.x
  211.             z=z+1
  212.             dat=data.s.x.y
  213.             sum.s.x=(dat-meen)**2+(sum.s.x)
  214.             if z=rows then do
  215.                 var.s.x=(sum.s.x)/N /* Variance & SD within cells */
  216.                 sd.s.x=sqrt(var.s.x)
  217.                 z=0
  218.                 s=s+1
  219.             end
  220.         end
  221. end
  222. summ.=0
  223. varrow.=0
  224. sdrow.=0
  225. s=1
  226. z=0
  227. N=((Ncols-1)*rows)-1
  228. do y=1 to NRows
  229.     z=z+1
  230.     do x=2 to NCols
  231.             meen=meanrow.s
  232.             dat=data.s.x.y
  233.             summ.s=(dat-meen)**2+(summ.s)
  234.     end
  235.     if z=rows then do
  236.         varrow.s=(summ.s)/N /* Variance & SD of sample */
  237.         sdrow.s=sqrt(varrow.s)
  238.         z=0
  239.         s=s+1
  240.     end
  241. end
  242. if fo then call writech(6Info,"40 ")
  243. summcol.=0
  244. varcol.=0
  245. sdcol.=0
  246. z=0
  247. N=NRows-1
  248. do x=2 to NCols
  249.     s=1
  250.     meen=meancol.x
  251.     do y=1 to NRows
  252.         z=z+1
  253.         dat=data.s.x.y
  254.         summcol.x=(dat-meen)**2+(summcol.x)
  255.             if z=rows then do
  256.                 z=0
  257.                 s=s+1
  258.             end
  259.     end
  260.     varcol.x=(summcol.x)/N /* Variance & SD of columns */
  261.     sdcol.x=sqrt(varcol.x)
  262. end
  263. /* Calculate Sum of Squares */
  264. SSR=0 /* Sum of Squares Rows */
  265. SSC=0 /* Sum of Squares Columns */
  266. SSW=0 /* Sum of Squares Within Cells */
  267. SSTotal=0 /* Sum of Squares Total */
  268. SSI=0 /* Sum of Squares Interaction */
  269.  
  270. /* Calculate SSR */
  271. ss=0
  272. do s=1 to NumSamples
  273.     ss=ss+((meanrow.s)-GM)**2
  274. end
  275. SSR=ss*(NCols-1)*rows
  276. /* Calculate SSC */
  277. ss=0
  278. do x=2 to NCols
  279.     ss=ss+((meancol.x)-GM)**2
  280. end
  281. SSC=ss*NumSamples*rows
  282. if fo then call writech(6Info,"50 ")
  283. /* Calculate SSI */
  284. ss=0
  285. do s=1 to NumSamples
  286.     do x=2 to NCols
  287.         ss=ss+((meanrc.s.x)-(meanrow.s)-(meancol.x)+GM)**2
  288.     end
  289. end
  290. SSI=ss*rows
  291.  
  292. /* Calculate SSW */
  293. ss=0
  294. z=0
  295. do x=2 to NCols
  296.     s=1
  297.         do y=1 to NRows
  298.             z=z+1
  299.             ss=ss+((data.s.x.y)-(meanrc.s.x))**2
  300.             if z=rows then do
  301.                 z=0
  302.                 s=s+1
  303.             end        
  304.         end
  305. end
  306. SSW=ss
  307. if fo then call writech(6Info,"60 ")
  308. /* Calculate SSTotal */
  309. ss=0
  310. z=0
  311. do x=2 to NCols
  312.     s=1
  313.         do y=1 to NRows
  314.             z=z+1
  315.             ss=ss+((data.s.x.y)-GM)**2
  316.             if z=rows then do
  317.                 z=0
  318.                 s=s+1
  319.             end
  320.         end
  321. end
  322. SSTotal=ss
  323. /* Calculate Variances */
  324. DFR=NumSamples-1 /* Degrees of Freedom Rows */
  325. DFC=NCols-2 /* Degrees of Freedom Columns */
  326. DFI=DFR*DFC /* Degress of Freedom Interaction */
  327. DFW=NumSamples*(NCols-1)*(Rows-1) /* Degrees of Freedom Within Cells */
  328. DFT=rows*NumSamples*(NCols-1)-1 /* Degrees of Freedom Total */
  329. VarR=SSR/DFR /* Variance Estimate for Samples */
  330. VarC=SSC/DFC /* Variance Estimate for Columns */
  331. VarW=SSW/DFW /* Variance estimate for Within Cells */
  332. VarI=SSI/DFI /* Variance Estimate for Interaction */
  333. Fcol=VarC/VarW /* F Ratio for Columns */
  334. Frow=VarR/VarW /* F Ratio for Rows */
  335. Finter=VarI/VarW /* F Ratio for Interaction */
  336. if fo then call writech(6Info,"70 ")
  337.  
  338. /* Calculate Probabilities */
  339.  
  340. AC=.0001
  341. EC=.005
  342. EC2=EC+EC
  343. PCol=PROBABILITY(Fcol,DFC,DFW)
  344. PCol=1-PCol
  345. PRow=PROBABILITY(Frow,DFR,DFW)
  346. PRow=1-PRow
  347. PInt=PROBABILITY(Finter,DFI,DFW)
  348. PInt=1-PInt
  349. /* P=.95 OR .99 */
  350. FCRIT1Col=FCRITICAL(.95,DFC,DFW)
  351. FCRIT2Col=FCRITICAL(.99,DFC,DFW)
  352. if fo then call writech(6Info,"80 ")
  353. /* P=.95 OR .99 */
  354. FCRIT1Row=FCRITICAL(.95,DFR,DFW)
  355. FCRIT2Row=FCRITICAL(.99,DFR,DFW)
  356. /* P=.95 OR .99 */
  357. if fo then call writech(6Info,"90 ")
  358. FCRIT1Int=FCRITICAL(.95,DFI,DFW)
  359. FCRIT2Int=FCRITICAL(.99,DFI,DFW)
  360. if fo then do
  361.     call writeln(6Info,"100 ")
  362.     call writeln(6Info,"Writing output to window...")
  363. end
  364. /* Output */
  365. 'SelectCell' out_cell
  366. 'ColumnWidth 10'
  367. 'Put "ANOVA - Two Way: With Replication"'
  368. 'CursorDown 2'
  369. 'Put "Summary Statistics"'
  370. 'CursorDown 2'
  371. 'GetCursorPos'
  372. go_cell=result
  373. do x=start_col+1 to end_col
  374.     'ColumnWidth 10'
  375.     'Column 1'
  376.     'Alignment 2'
  377.     title=""""||title.x||""""
  378.     'Put' title
  379. end
  380. 'Column 1'
  381. 'ColumnWidth 10'
  382. 'Alignment 2'
  383. 'Put' "Total"
  384. 'SelectCell' go_cell
  385. 'CursorDown 2'
  386. do s=1 to NumSamples
  387.     'GetCursorPos'
  388.     start_cell=result
  389.     sampletitle=""""||sampletitle.s||""""
  390.     'Put' sampletitle
  391.     'CursorDown 1'
  392.     'Put' "Count:"
  393.     'CursorDown 1'
  394.     'Put' "Sum:"
  395.     'CursorDown 1'
  396.     'Put' "Mean:"
  397.     'CursorDown 1'
  398.     'Put' "Variance:"
  399.     'CursorDown 1'
  400.     'GetCursorPos'
  401.     end_cell=result
  402.     'SelectCell' start_cell
  403.     'CursorDown 1'
  404.     do x=2 to NCols
  405.         'Column 1'
  406.         'Put' rows
  407.         'CursorDown 1'
  408.         'Put' cell.s.x
  409.         'CursorDown 1'
  410.         'Put' format(meanrc.s.x,,4)
  411.         'CursorDown 1'
  412.         'Put' format(var.s.x,,4)
  413.         'CursorUp 3'
  414.     end
  415.     'Column 1'
  416.     'Put' scount
  417.     'CursorDown 1'
  418.     'Put' rowtot.s
  419.     'CursorDown 1'
  420.     'Put' format(meanrow.s,,4)
  421.     'CursorDown 1'
  422.     'Put' format(varrow.s,,4)
  423.     'CursorDown 1'
  424.     'SelectCell' end_cell
  425.     'CursorDown 2'
  426. end
  427. 'Put "Col. Total"'
  428. 'CursorDown 2'
  429. 'GetCursorPos'
  430. start_cell=result
  431. 'Put' "Count"
  432. 'CursorDown 1'
  433. 'Put' "Sum:"
  434. 'CursorDown 1'
  435. 'Put' "Mean:"
  436. 'CursorDown 1'
  437. 'Put' "Variance:"
  438. 'GetCursorPos'
  439. end_cell=result
  440. 'SelectCell' start_cell
  441. do x=2 to NCols
  442.     'Column 1'
  443.     'Put' count.x
  444.     'CursorDown 1'
  445.     'Put' total.x
  446.     'CursorDown 1'
  447.     'Put' format(meancol.x,,4)
  448.     'CursorDown 1'
  449.     'Put' format(varcol.x,,4)
  450.     'CursorUp 3'
  451. end
  452. 'SelectCell' end_cell
  453. 'CursorDown 2'
  454. 'Put' "ANOVA"
  455. 'CursorDown 2'
  456. 'GetCursorPos'
  457. start_cell=result
  458. 'Alignment 2'
  459. 'Put "Source of"'
  460. 'CursorDown 1'
  461. 'Alignment 2'
  462. 'Put "Variation"'    
  463. 'Column 1'
  464. 'CursorUp 1'
  465. 'Alignment 2'
  466. 'Put "Sum of"'
  467. 'CursorDown 1'
  468. 'Alignment 2'
  469. 'Put "Squares"'
  470. 'Column 1'
  471. 'CursorUp 1'
  472. 'Alignment 2'
  473. 'Put "Degrees of"'
  474. 'CursorDown 1'
  475. 'Alignment 2'
  476. 'Put "Freedom"'
  477. 'Column 1'
  478. 'CursorUp 1'
  479. 'Alignment 2'
  480. 'Put' "Variance"
  481. 'CursorDown 1'
  482. 'Alignment 2'
  483. 'Put "Estimate"'
  484. 'Column 1'
  485. 'CursorUp 1'
  486. 'Alignment 2'
  487. 'Put "F Ratio:"'
  488. 'SelectCell' Start_cell
  489. 'CursorDown 2'
  490. 'Put' "Sample:"
  491. 'CursorDown 1'
  492. 'Put' "Column:"
  493. 'CursorDown 1'
  494. 'Put' "Interaction:"
  495. 'CursorDown 1'
  496. 'Put "Within Cells:"'
  497. 'CursorDown 1'
  498. 'Put' "Total:"
  499. 'CursorDown 2'
  500. 'GetCursorPos'
  501. test_cell=result
  502. 'SelectCell' Start_cell
  503. 'CursorDown 2'
  504. 'Column 1'
  505. 'Put' format(SSR,,4)
  506. 'CursorDown 1'
  507. 'Put' format(SSC,,4)
  508. 'CursorDown 1'
  509. 'Put' format(SSI,,4)
  510. 'CursorDown 1'
  511. 'Put' format(SSW,,4)
  512. 'CursorDown 1'
  513. 'Put' format(SSTotal,,4)
  514. 'CursorUp 4'
  515. 'Column 1'
  516. 'Put' DFR
  517. 'CursorDown 1'
  518. 'Put' DFC
  519. 'CursorDown 1'
  520. 'Put' DFI
  521. 'CursorDown 1'
  522. 'Put' DFW
  523. 'CursorDown 1'
  524. 'Put' DFT
  525. 'CursorUp 4'
  526. 'Column 1'
  527. 'Put' format(VarR,,4)
  528. 'CursorDown 1'
  529. 'Put' format(VarC,,4)
  530. 'CursorDown 1'
  531. 'Put' format(VarI,,4)
  532. 'CursorDown 1'
  533. 'Put' format(VarW,,4)
  534. 'CursorUp 3'
  535. 'Column 1'
  536. 'Put' format(Frow,,4)
  537. 'CursorDown 1'
  538. 'Put' format(Fcol,,4)
  539. 'CursorDown 1'
  540. 'Put' format(Finter,,4)
  541. 'SelectCell' test_cell
  542. 'Put "P(F Sample <=f):"'
  543. 'CursorDown 1'
  544. 'Put "F-Critical One-Tail(95%):"'
  545. 'CursorDown 1'
  546. 'Put "F-Critical One-Tail(99%):"'
  547. 'CursorDown 1'
  548. 'Put "P(F Column <=f):"'
  549. 'CursorDown 1'
  550. 'Put "F-Critical One-Tail(95%):"'
  551. 'CursorDown 1'
  552. 'Put "F-Critical One-Tail(99%):"'
  553. 'CursorDown 1'
  554. 'Put "P(F Interaction <=f):"'
  555. 'CursorDown 1'
  556. 'Put "F-Critical One-Tail(95%):"'
  557. 'CursorDown 1'
  558. 'Put "F-Critical One-Tail(99%):"'
  559. 'SelectCell' test_cell
  560. 'Column 2'
  561. 'Put' format(PRow,,6)
  562. 'CursorDown 1'
  563. 'Put' format(FCRIT1Row,,4)
  564. 'CursorDown 1'
  565. 'Put' format(FCRIT2Row,,4)
  566. 'CursorDown 1'
  567. 'Put' format(PCol,,6)
  568. 'CursorDown 1'
  569. 'Put' format(FCRIT1Col,,4)
  570. 'CursorDown 1'
  571. 'Put' format(FCRIT2Col,,4)
  572. 'CursorDown 1'
  573. 'Put' format(PInt,,6)
  574. 'CursorDown 1'
  575. 'Put' format(FCRIT1Int,,4)
  576. 'CursorDown 1'
  577. 'Put' format(FCRIT2Int,,4)
  578. 'Refresh 1'
  579. 'Refresh 2'
  580. /*'Message' "Finished"*/
  581. /*indicate the main script is finished*/
  582. DisplayMsg="Cleaning up ...."||CR||"Exiting"
  583. result=writeln(6Info, DisplayMsg)
  584. if result~=0 then do
  585.     /*Wait 3 seconds*/
  586.     CALL DELAY(150)
  587.     /* close window*/
  588.     result=close(6Info)
  589. end
  590. 'DEFPUBSCREEN("Workbench")'
  591. exit
  592.  
  593. /* Procedures */
  594.  
  595. cellrow: procedure
  596. do
  597.     parse arg cell
  598.     do charpos=2 to length(cell)
  599.     if datatype(substr(cell,charpos,1),n) then return substr(cell,charpos)
  600.     end
  601.     return 0
  602. end
  603. Return
  604.  
  605. cellcol: procedure
  606. do
  607.     parse arg cell
  608.     labels="ABCDEFGHIJKLMNOPQRSTUVWXYZ"
  609.     cell=upper(cell)
  610.     len=length(cell)
  611.     val=0
  612. do charpos=1 to len
  613.     if datatype(substr(cell,charpos,1),n) then
  614.     do cell=reverse(substr(cell,1,charpos-1))
  615.     do x=1 to length(cell)
  616.     val=(26**(x-1))*pos(substr(cell,x,1),labels)+val
  617.     end
  618.     return val
  619.     end
  620.     end
  621.     return 0
  622. end
  623. Return
  624.  
  625. syntax:
  626.      if arg(1)='FAIL' then do
  627.         'Message' "Library is unavailable."
  628.         'DEFPUBSCREEN("Workbench")'
  629.         exit
  630.         end
  631.     'Message' "Unknown Syntax Error"
  632.     'DEFPUBSCREEN("Workbench")'
  633.     exit
  634.  
  635. Format:  procedure
  636.  
  637.      arg number, before, after
  638.      CallLine = SIGL
  639.      if ~datatype(CallLine, 'N') then CallLine = '??'
  640.  
  641.      /* Make sure we have a number as first (required) argument    */
  642.      if ~datatype(number, 'N') then do
  643.         if number = '' then
  644.            fc = 17     /* Wrong number of arguments           */
  645.         else
  646.            fc = 47     /* Arithmetic conversion error             */
  647.         signal FormatSyntaxError
  648.      end
  649.      num = number + 0
  650.      if before = '' & after = '' then
  651.         return num
  652.      else do
  653.         parse var num integer '.' fraction
  654.         if before = '' then before = length(integer)
  655.         if after = '' then after = length(fraction)
  656.         if ~datatype(before, N) | ~datatype(after, N) then
  657.            do fc = 18
  658.            signal FormatSyntaxError
  659.        end
  660.         if before < length(integer) then do
  661.            fc = 18
  662.            signal FormatSyntaxError
  663.         end
  664.         if after ~= length(fraction) then do
  665.            fraction = trunc(('.'fraction'0') + ('.'copies('0', after)'5'), after)
  666.         if integer<1&integer>-1 then integer=integer
  667.            else integer = integer + (fraction % 1)
  668.            fraction = substr(fraction, 3)
  669.         end
  670.         if fraction >= 0 then
  671.            return right(integer, before)'.'fraction
  672.         else
  673.            return right(integer, before)
  674.      end
  675.  
  676.  FormatSyntaxError:
  677.         if show('F', STDERR) then
  678.            call writeln(STDERR, '+++ Error' fc 'in line' CallLine':' errortext(fc))
  679.         else
  680.            mess='+++ Error' fc 'in line' CallLine':' errortext(fc)
  681.         'Message' mess
  682.         parse source Func .
  683.         if Func = 'FUNCTION' then do
  684.         'DEFPUBSCREEN("Workbench")'
  685.            exit "Err"
  686.         end
  687.         else do
  688.         'DEFPUBSCREEN("Workbench")'
  689.            exit 10
  690.         end
  691.  
  692. FCRITICAL: PROCEDURE EXPOSE AC EC EC2
  693.  
  694.     PARSE ARG PO,K1,K2
  695.     /* fIND NORMAL Z CORRESPONDING TO P */
  696.     P1=PO
  697.     Z=NORMALPP(PO)
  698.     IF K2<3 THEN DO
  699.         W=K2**.75
  700.         W=Z/W
  701.         Z=Z*(1.1581-W*(.2296+W*(.0042+.0027*W)))
  702.     END
  703.     /* FIND INITIAL APPROX. TO F */
  704.     A1=2/(9*K1)
  705.     A2=2/(9*K2)
  706.     W=Z*Z
  707.     W1=1+A2*(A2-W-2)
  708.     W2=A1+A2-A1*A2-1
  709.     W3=1+A1*(A1-W-2)
  710.     SN=SIGN(W2*W2-W1*W3)
  711.     IF SN=-1 THEN W3=-1*SQRT(ABS(W2*W2-W1*W3))
  712.     ELSE W3=SQRT(W2*W2-W1*W3)
  713.     FO=(W3-W2)/W1
  714.     FO=FO**3
  715.     /* MODIFIED NEWTON ITERATION TO IMPROVE VALUE OF F */
  716.     DO FOREVER
  717.         FCRIT=FO
  718.         PO=PROBABILITY(FCRIT,K1,K2)
  719.         F1=PO-P1
  720.         FCRIT=FO+EC
  721.         PO=PROBABILITY(FCRIT,K1,K2)
  722.         F2=PO
  723.         FCRIT=FO-EC
  724.         PO=PROBABILITY(FCRIT,K1,K2)
  725.         F2=F2-PO
  726.         F2=F2/EC2
  727.         F1=FO-F1/F2
  728.         IF ABS(F1-FO)>AC THEN FO=F1
  729.         ELSE LEAVE
  730.     END
  731.     FCRIT=F1
  732. RETURN FCRIT
  733.  
  734. NORMALPP: PROCEDURE
  735.  
  736.     ARG P0
  737.     A1=2.515517
  738.     A2=.802853
  739.     A3=.010328
  740.     A4=1.432788
  741.     A5=.189269
  742.     A6=.001308
  743.     Q=.5-ABS(P0-.5)
  744.     SN=SIGN(-2*LN(Q))
  745.     IF SN=-1 THEN W=-1*SQRT(ABS(-2*LN(Q))
  746.     ELSE W=SQRT(-2*LN(Q))
  747.     W1=A1+W*(A2+A3*W)
  748.     W2=1+W*(A4+W*(A5+A6*W))
  749.     ZZ=W-W1/W2
  750.     SELECT
  751.         WHEN (P0-.5)<0 THEN TT=-1
  752.         WHEN (P0-.5)=0 THEN TT=0
  753.         otherwise TT=1
  754.     END
  755.     ZZ=ZZ*TT
  756. RETURN ZZ
  757.  
  758. PROBABILITY: PROCEDURE EXPOSE AC EC EC2
  759.  
  760.     PARSE ARG F,K1,K2
  761.     IF K1=1 THEN AC=.00001
  762.     H2=K2/2
  763.     H1=K1/2
  764.     H3=H1+H2
  765.     L1=0
  766.     XX=H1
  767.     LG=LOGGAMMA(XX)
  768.     L1=L1+LG
  769.     XX=H2
  770.     LG=LOGGAMMA(XX)
  771.     L1=L1+LG
  772.     XX=H3
  773.     LG=LOGGAMMA(XX)
  774.     L1=L1-LG
  775.     W1=K2/(K2+K1*F)
  776.     XX=0
  777.     IF H2<(H3*W1) THEN DO
  778.         W2=W1
  779.         W1=1-W1
  780.         W3=H2
  781.         H2=H1
  782.         H1=W3
  783.     END
  784.     ELSE DO
  785.         W2=1-W1
  786.         XX=1
  787.     END
  788.     T1=1
  789.     W3=1
  790.     P=1
  791.     I=INT(H1+W2*H3)
  792.     W4=W1/W2
  793.     /*LABELB:*/
  794.     T2=H1-W3
  795.     IF I=0 THEN W4=W1
  796.     /*LABELC:*/
  797.     DO FOREVER
  798.         T1=T1*T2*W4/(H2+W3)
  799.         P=P+T1
  800.         T2=ABS(T1)
  801.         IF (T2<=AC) & (T2<=AC*P) THEN LEAVE /*SIGNAL 'LABELD'*/
  802.         W3=W3+1
  803.         I=I-1
  804.         IF I>=0 THEN DO /*SIGNAL 'LABELB'*/
  805.             T2=H1-W3
  806.             IF I=0 THEN W4=W1
  807.             ITERATE
  808.         END
  809.         T2=H3
  810.         H3=H3+1
  811.     /*SIGNAL 'LABELC'*/
  812.     END
  813.     /*LABELD:*/
  814.     W1=EXP(H2*LN(W1)+(H1-1)*LN(W2)-L1)
  815.     P=P*W1/H2
  816.     IF XX=0 THEN PO=P
  817.     ELSE PO=1-P
  818.     AC=.001
  819. RETURN PO
  820.  
  821. LOGGAMMA: PROCEDURE
  822.  
  823.     ARG XA
  824.     C1=76.18009173
  825.     C2=-86.50532033
  826.     C3=24.01409822
  827.     C4=-1.231739516
  828.     C5=.001208580
  829.     C6=-.000005364
  830.     C7=2.506628275
  831.     X1=XA-1
  832.     W=X1+5.5
  833.     W=(X1+.5)*LN(W)-W
  834.     S=1+C1/(X1+1)+C2/(X1+2)+C3/(X1+3)
  835.     S=S+C4/(X1+4)+C5/(X1+5)+C6/(X1+6)
  836.     L=W+LN(C7*S)
  837. RETURN L